function [DFILE, DPARS, err] = ana_setup(DFILE, sf)
% general setup for analysis...
% 6/2/2002 P. Manis
% modified 7/19/2007 P. Manis
% now also handles the psp times etc...
%


global VOLTAGE CURRENT ALLCH
global CONTROL

DPARS = [];
err = 1;
%-----------------------------General Setup--------------------------------

[DFILE1] = synch_file(DFILE, sf);
% do the "block info" thing to get stimulus parameters loaded into dfile.
%[sf, df] = block_info(DFILE1.fullfile, DFILE1.dblock);


if(isempty(DFILE1) | isempty(VOLTAGE))
    QueMessage('Analysis Setup: could not get DFILE');
    return;
end;
DFILE = DFILE1;

% verify fields in control structure. Create if they don't exist and then return.
update = 0; % flag
if(~check_field(CONTROL(sf),'stim_time')) % doesn't exist - then create it.
    CONTROL(sf).stim_time=[5 15 25 35 45];
    %	disp('stim_time')
    update = 1;
end;
if(~check_field(CONTROL(sf),'psp_time'))
    CONTROL(sf).psp_time= [5.5 14 15.5 24 25.5 34 35.5 44 45.5 54];
    %	disp('psp_time')
    update=1;
end
%if(~check_field(CONTROL(sf),'psp_reclist'))
%	CONTROL(sf).psp_reclist= ' ';
%	disp('reclist')
%	update=1;
%end

% check for units and insert them if they are not present.
if(~check_field(CONTROL(sf), 'I_Unit') | isempty(CONTROL(sf).I_Unit))
    CONTROL(sf).I_Unit='pA';
end

if(~check_field(CONTROL(sf), 'V_Unit') | isempty(CONTROL(sf).V_Unit))
    CONTROL(sf).V_Unit='mV';
end

if(~check_field(CONTROL(sf), 'T_Unit') | isempty(CONTROL(sf).T_Unit))
    CONTROL(sf).T_Unit='ms';
end
if(~isnumeric(CONTROL(sf).vgain)) | isempty(CONTROL(sf).vgain)
    CONTROL(sf).vgain = str2num(CONTROL(sf).vgain);
end;
if(~isnumeric(CONTROL(sf).igain) | isempty(CONTROL(sf).igain))
    CONTROL(sf).igain = str2num(CONTROL(sf).igain);
end;

if(~check_field(CONTROL(sf), 'sign_list') | isempty(CONTROL(sf).sign_list))
    CONTROL(sf).sign_list=[-1 -1];
end;


DPARS = ana_setup2(DFILE, sf); % get basic information and record times

% get the step windows
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);
tsw3 = 0;
tsw4 = 0;
hypflag = 0;
all_ivs = {'cciv', 'cciv-l', 'ap-iv', 'ap_iv', 'ap-iv2', 'ap-ivl', 'ap-ivl2', 'vc_iv', 'vciv', 'vc_act'};
all_hyps = {'hyp', 'hyps', 'hypsave', 'hypl', 'hypwatch', 'ap-hyp', 'ap_hyp', ...
    'ap_hypsynch', 'ap-hyp2', 'ap-hyp3', 'mc-hyp2', 'ltp_base'};
DPARS.tdel = ts1;
DPARS.tstepdur = ts2;

% designate times according to protocol
switch(DPARS.protocol)
    case all_ivs % standard IV from -60 mV
        tdel=ts1; tstepdur=ts2;
        tsl = ts1; % spike latency reference point
        tsw1 = ts1;
        tsw2 = ts2;
    case 'iv'
        tdel=ts1; tstepdur=ts2;
        tsl = ts1;
        tsw1 = ts1;
        tsw2 = ts2;
    case 'ivf'
        tdel=ts1; tstepdur=ts2;
        tsl = ts1;
        tsw1 = ts1;
        tsw2 = ts2;
    case 'drugmonitor'
        tdel = DFILE.durs(1);
        tstepdur = DFILE.durs(4);
        ts1 = tdel + DFILE.durs(2);
        tsw1=ts1;
        tsw2 = tsw1 + DFILE.durs(3);
        tsw3 = tsw2 + DFILE.durs(4);
        tsw4 = tsw3 + DFILE.durs(5);


    case all_hyps % variable prepulse level before evoking spikes
        tstepdur=ts2; % measure IV during prepulse
        tsw1 = ts1+ts2;
        tdel = ts1;
        tsw2 = ts1+ts2+ts3; % but do spike analysis/counts in second window
        tsl = tsw1;
        hypflag = 1; % set so we can do boltzman fit...
    otherwise % just dump some defaults, - although this might cause problems later on
        DPARS.tdel = ts1;
        DPARS.tsl = ts1+ts2;
        DPARS.tsw1 = ts1;
        DPARS.tsw2 = ts2;
        DPARS.tsw3 = ts3;
        DPARS.tsw4 = ts4;
        DPARS.tstepdur = ts2;
        DPARS.hypflag = hypflag;
        %        QueMessage(sprintf('ana_setup - Error: protocol %s not recoginized - Exiting', DPARS.protocol));
        err = 0;

end

tsa = [ts1, ts1+ts2, ts1+ts2+ts3, ts1+ts2+ts3+ts4];

DPARS.tsw1 = tsa(1); % thse are actual times - use tstep for indices.
DPARS.tsw2 = tsa(2);
DPARS.tsw3 = tsa(3);
DPARS.tsw4 = tsa(4);
DPARS.hypflag = hypflag;
DPARS.tsw = tsa; % save the array of times as well.

% calculate point indices for each step (ho, s1, s2, s3).
for i = 1:length(tsa)
    DPARS.tstep(i,:) = floor(tsa(i)./DPARS.RATES-0.5); % pointer to end of durho
    DPARS.tstep(find(DPARS.tstep > DPARS.pts)) = DPARS.pts;  % make sure steps are always in a valid range
    DPARS.tstep(find(DPARS.tstep < 1)) = 1;  % make sure steps are always in a valid range 
end;

% calculate times and point indices for psp windows
QueMessage('Setting PSP/PSC analysis windows');
[DPARS.tpsp, DPARS.tpspi, DPARS.psp_pair] = parse_array(DPARS, CONTROL(sf).psp_time);
[DPARS.tstim, DPARS.tstimi, DPARS.stim_pair] = parse_array(DPARS, CONTROL(sf).stim_time);

QueMessage('ana_setup successful & complete', 1);
err = 0;
return;




function [ts, tp, pair] = parse_array(DPARS, atime)
% parses psp and stim arrays from teh control list
% returns the times (in ts) and the pointers (in tp)
% as 2d and 3d arrays.
% ts contains the times in an nxm array, just as the data is
% tp contains the points in an nxmxk array, wehre each k corresponds to a
% successive record (allowing for records to have different sampling times
%
if(ischar(atime))
    atime = str2num(atime); % convert to numbers if it is a string
end;
l = size(atime);
if(prod(l) == 4)
     % reshape for multi stimulus...
    atime = reshape(atime, 1, 4)'; % multistimulus versions
end
if(prod(l) == 8)
    atime = reshape(atime, 1, prod(l))';
end;
p2=size(atime);
pair=(p2(1)>1);

[nipsp, nkpsp] = size(atime);

ts = zeros(nipsp, nkpsp);
for i = 1:nipsp
    for k = 1:nkpsp
        ts(i,k)=number_arg(atime(i,k)); % get the baseline dur
    end;
end;

% calculate measurement times

for i = 1:nipsp
    for k = 1:nkpsp
        tp(i,k,:) = floor(ts(i,k)./DPARS.RATES - 0.5); % for RMP/ihold determination
    end;
end;





%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

function [DPARS] = ana_setup2(DFILE, sf)
% setup routine to extract the time information:
% gets the time array for each data record (computed)
% and gets the sample
global CONTROL ALLCH

time = [];

% abstract general information
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;
[records,pts]=size(ALLCH{1});

[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

%QueMessage('Analysis setup:  time base');
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE)';
tmax=max(max(time));
%k = find(DFILE.rate < 1);
%DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec

% always compute the time of each stimulus, in seconds from the start of the data
%QueMessage('Analysis setup:  time base 2');
if(DFILE.mode >= 5)
    wz=DFILE.ztime; %(RL-RL(1)+1); % indexing allows for skip parameters... etc.
    w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
    if(~isempty(w))
        wz(w+1:end)=wz(w+1:end)+24*60*60;
    end;
    zt = (wz-wz(1))/(60);
    cond_baseline = 5; % 5 min baseline
else
    zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
    cond_baseline = 3;
end;
%TM=zt;
TM=DFILE.ztime/60; %(RL-RL(1)+1);

%QueMessage('Analysis setup:  time base 3');
if(DFILE.mode >= 5)
    wz=DFILE.ztime;
    w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
    if(~isempty(w))
        wz(w+1:end)=wz(w+1:end)+24*60*60;
    end;
    ZT = (wz-wz(1))/(60);
    cond_baseline = 5; % 5 min baseline
elseif length(DFILE.ztime) == length(RL)
    ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
    cond_baseline = 3;
else
    ZT = ones(length(RL), 1)*DFILE.cycle;
    cond_baseline = 5;
end;
uz = find(diff(ZT) < 0); % look for time reversals - counter was reset. If found, add last time + stim interval and restart.
if(~isempty(uz))
    lj = 1;
    for j = uz
        ZT(j+1:end) = ZT(j+1:end) + ZT(j)+(ZT(j)-ZT(j-1));
    end;
end;

% access ZT with the record number (1..n) to get the corresponding ztime
% have to get the records first

% Now, get the times when the valves switched (if any...)
% and generate periods with valve 1, 2, 3 or 4.
%QueMessage('Analysis setup:  Valves...');
p=datac('getnote'); % read the current notefile information.
t_sw_valve=[];
sw_valve=[];
n_valve=[];
TL=[];
VL=[];
if(~isempty(p) && ~isempty([p.proto])) % there should be some, but if not, don;'t do much
    % first set of arrays are immediate representations.
    sw_valve=[1 [p(find(diff([p.valve])~=0)+1).frec]];	% valve switch list (records)
    if(length(sw_valve) > 1)
        n_valve=[1 p(find(diff([p.valve])~=0)+1).valve]; % which valve...
        t_sw_valve=ZT(sw_valve); % don't forget offset from start of data..
        % now make long time arrays to match the other arrays.
        for i = frec:lrec
            TL(i)=ZT(i);
            for j=1:length(sw_valve)
                if(i >= sw_valve(j))
                    VL(i)=n_valve(j);
                end;
            end;
        end;
    end;
end;
DPARS.protocol = protocol;
DPARS.ftime = DFILE.ftime;
DPARS.rate = rate;
DPARS.RATES = RATES;
DPARS.records = records;
DPARS.pts = pts;
DPARS.frec = frec;
DPARS.lrec = lrec;
DPARS.spike_thresh=number_arg(CONTROL(sf).thresh);
DPARS.time = time';
DPARS.TM = TM;
DPARS.TL = TL;
DPARS.ZT = ZT;
DPARS.VL = VL;
DPARS.cond_baseline = cond_baseline;
DPARS.n_valve = n_valve;
DPARS.sw_valve = sw_valve;
